rm(list = ls())
gc()
require(tidyverse)
# require(marginaleffects)
# require(patchwork)
require(mcp)

################################################################################
##
## Date:    2024-12-22
## Author:  james.h.bisbee@vanderbilt.edu
## Purpose: This script generates SI Figure 7.
## Inputs:  /scratch/jhb362/zilinsky_2023/data/results/VIMP_ranger/VIMP_2024_outcome-.*(CUTBACKSPEND|MON|FEEL|SPEND).*-weeks.*_DumFacts-FALSE_temp-chg12.RData
##            - Variable importance results generated by NFR_vimp_prep.R
##            - Summarized on the NYU HPC into PSRM_comb_weeks_chg12.RData via NFR_data_prep.R
## Outputs: ./figures/figS7.pdf
## Notes:   For this to run, JAGS version 4 is required.
##
################################################################################

# Compute details
print(paste0('Compute environment from ',Sys.Date(),' run by Bisbee'))
if(Sys.info()['sysname'] == 'Windows') {
  ram_size = system("wmic MemoryChip get Capacity", intern = TRUE)[-1]
  model_name = system("wmic cpu get name", intern = TRUE)[2] # nocov
  vendor_id = system("wmic cpu get manufacturer", intern = TRUE)[2] # nocov
  
  print(list(ram = stringr::str_squish(ram_size)[1],
             vendor_id = stringr::str_squish(vendor_id),
             model_name = stringr::str_squish(model_name),
             no_of_cores = parallel::detectCores()))
} else if(Sys.info()['sysname'] == 'Linuxs') {
  splitted <- strsplit(system("ps -C rsession -o %cpu,%mem,pid,cmd", intern = TRUE), " ")
  df <- do.call(rbind, lapply(splitted[-1], 
                              function(x) data.frame(
                                cpu = as.numeric(x[2]),
                                mem = as.numeric(x[4]),
                                pid = as.numeric(x[5]),
                                cmd = paste(x[-c(1:5)], collapse = " "))))
  df
} else {
  cat("If not on Linux or Windows, you'll have to figure out your own solution to seeing the compute environment.")
}

sessionInfo()

load('./data/VIMP_ranger/PSRM_comb_weeks_chg12.RData')

toplot <- vimp %>%
  mutate(relVimp = vimp / predErr) %>%
  select(vars,outcome,period,n,relVimp) %>%
  group_by(vars,outcome,period) %>%
  summarise_all(list(mean = ~mean(.,na.rm=T),
                     lb = ~quantile(.,.005,na.rm=T),
                     ub = ~quantile(.,.995,na.rm=T))) %>%
  ungroup() %>%
  left_join(lookup)

model = list(
  relVimp_mean ~ 1,
  ~ 1 + period2
)



plts <- list()
for(y in toplot %>% count(outcome) %>% pull(outcome)) {
  # stop()
  fit = mcp(model,data = toplot %>%
              mutate(period = as.Date(period)) %>%
              filter(period >= as.Date('2010-01-01'),
                     period <= as.Date('2016-01-01')) %>%
              filter(vars == 'PARTY',
                     outcome == y) %>%
              arrange(period) %>%
              mutate(period2 = as.numeric(period)) %>%
              data.frame(),
            par_x = "period2")
  
  plts[[y]] <- plot(fit) + 
    scale_x_continuous(labels = function(x) as.Date(x,origin = '1970-01-01'),
                       breaks = as.numeric(as.Date(c('2012-01-01','2012-07-01','2013-01-01','2013-07-01')))) + 
    geom_vline(xintercept = as.numeric(as.Date('2013-01-01')),
               linetype = 'dashed')
}


pdf('./figures/figS7.pdf',width = 8,height = 8)
((plts$CUTBACKSPEND & 
    labs(x = NULL,y = NULL,
         title = "Party ID Variable Importance by Week",
         subtitle = 'Cutting back on spending') & 
    theme_bw()) + (plts$FEELBETTERFIN & 
                     labs(x = NULL,y = NULL,
                          subtitle = 'Feeling better about finances') & 
                     theme_bw())) / 
  ((plts$FEELGOODMON & 
      labs(x = NULL,y = 'Party ID Variable Importance',
           subtitle = 'Feeling good about money') & 
      theme_bw()) + (plts$MONMORETHANENOUGH & 
                       labs(x = NULL,y = NULL,
                            subtitle = 'More than enough money') & 
                       theme_bw())) / 
  ((plts$MONWORRYYEST & 
      labs(x = 'Date',y = NULL,
           subtitle = 'Worried about money') & 
      theme_bw()) + (plts$WATCHSPEND & 
                       labs(x = 'Date',y = NULL,
                            subtitle = 'Watching spending') & 
                       theme_bw()))
dev.off()

# EOF